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ABSTRACT 

The mass function of dark halos in a A-dominated cold dark matter (ACDM) universe is investigated. 529 
output files from five runs of A^-body simulations are analyzed using the friends-of-friends cluster finding 
algorithm. All the runs use 512 3 particles in the box size of 35 /r'Mpc to 140 /r'Mpc. Mass of particles for 
35 h~ l Mpc runs is 2.67 x 10 7 /z~' Mq. Because of the high mass resolution of our simulations, the multiplicity 
function in the low-mass range, where the mass is well below the characteristic mass and v = S c /a < 1.0, is 
evaluated in the present work, and is well fitted by the functional form proposed by Sheth & Tormen (ST). 
However, the maximum value of the multiplicity function from our simulations at v ~ 1 is smaller, and its low 
mass tail is shallower when compared with the ST multiplicity function. 

Subject headings: cosmology: theory — dark matter — galaxies: luminosity function, mass function 



1. INTRODUCTION 

The mass function is one o f the most import a nt st atisti- 
cal quantities of dark halos. iPress & SchechteJ Jl974ft in- 
vented an ansatz to predict the mass function of dark ha- 
los formed through hierarchical clustering, assuming that 
each smoothed mass element with arbitrary smoothing length 
evolves independently, in accordance with the spherical 
top-hat model. The PS formalism was exten ded so that 
the merging history of each halo is trac eable i lBond et alJ 
(19911 iBowetl [T99"U ILacev & Cold" Tl 993). This extended 
PS formalism has been used in a wide range of applica- 
tio ns. Especially the so-called semi-analytic galaxy mod- 



els (Cole et al. 2000; Kauffmann, White, & Guiderdoni 1992 , 
Nagashima et alJl200lU Somerville & Primack 1999) use it to 
successfully reproduce many properties of observed galaxies 
The PS mass function and the mass function from N- 
body s imulatio ns agre e with each other only qualitatively 
(Brain erd & Villumsenl fl99l lEfst athiou et al. 1988), and 
modification to threshold linear overdensity leads to a bet- 
ter agreement (Efstathiou & Rees 1988; Gelb & Bertschinger 
119941 ILacev & Cold 119941) . However, the PS mass 
function giv es the smaller number o f high -mass halos 
( iGross et alJ 119981 Uain & BertschingeJ 11994 . while giv- 
ing the larger number of low-mass hal os when com- 
pared with the jV-body mas s functions (iGovernato et all 
1999; Somervi lle et all 12000). This tendency was con- 
firmed over a mass range by connecting the numerical 
mass functions from simulations wit h different box sizes 
llGross et alJll998t iTenkins et alJl2001l) . Taking i nto account 
the s p atial correlation of density fluctu ations (Nagashima 
2001] lYano. N agashima. & Gouda 1996), or incorporating 
the ellipsoidal collapse model into the PS ansatz instead 
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of the spherical collapse model (lAu diLTevssier^&AJimj 
1 19971 lEnsteinlll984l ITee & Shandarin 1998; Monaco i}l995t 
Sheth. Mo! & Tormenl 120011 ISheth & Tormenl 120021) . an 
agreement was improved between the numerical mass func- 
tion and the analytic mass function, especially the ST mass 
function proposed by Sheth & Tormen ( 1999). In addition, 
the mass function of clus ter progenitors has also been studied 
(Okarnoto&Haberf999l). 

On the other hand, the numerical mass function depends on 
the c lus ter finding algori th m ado pt ed (iGelb & Bert schingerl 
1 19941 IGovernato et alJ 119991 ILacev & Cold 119931) . 
Jenki nset al.1 d20Qll hereafter J01) demonstrated that 

( iDavis et all 



the friends-of-friends (FoF) algorithm ( IDavis et al.1 1985 
with a fixed linking length of 0.2 times the mean particle sep- 
aration results in the numerical mass function that shows the 
best universality for various compositions of cosmological 
parameters and box sizes. Thus, it is better to set the threshold 
density proportional to the b ackground ma ss density rather 
than to the critical density. IWhitd d2002l hereafter W02) 
supported the above argument and demonstrated that an 
other definition of mass of halos, M^Qb, gives the universal 
numerical mass function. Here, M\^b is the mass within a 
sphere whose average density is 180 times the background 
density. 

However, different authors support different mass func- 
tions. For example, J01 propos ed a ne w fitting formula to the 
universal mass fu nction, while | W02l supported the ST mass 
function (see also Reed et al. 2003). Since the mass function 
of high-mass halos has extensively been studied by them, this 
discrepancy is possibly due to the lack of information on the 
behavior of the numerical mass function for low-mass halos. 
Thus, in order to investigate the functional form of the univer- 
sal mass function, we performed five runs of A^-body simula- 
tions with high mass resolution. 

In this paper, we briefly describe our simulation code and 
cosmological and simulation parameters adopted in §2. Re- 
sults of the simulations as well as the parameter values for 
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the best-fit mass function are given in §3, and discussions are 
given in §4. 

2. MULTIPLICITY FUNCTION 

The multiplicity function is the differential distribution 
function of the normalized fluctuation amplitude of dark ha- 
los for each mass elem ent. According to the definition by 
ISheth & Tormenl 11999). the multiplicity function is defined 

as 



vf(v) = M 



2 n(M,z) dlogM 



d\ogv 



(1) 



where n(M,z) is the number density of dark halos, v = 5 c /<jm 
is the peak height of a halo, 8 C is the linear overdensity at 
the collapse epoch of halos given by the spherical collapse 
model, and ctm is the standard deviation of the density fluc- 
tuation field smoothed by the top-hat window function. Note 
that the time evolution of the number density and the depen- 
dence of the number density on the initial power spectrum 
are absorbed in v. Since this universality of the multiplicity 
function is guaranteed under the PS ansatz, we compare our 
simulation results with analytic multipl icity functions. 

Along with the PS ansatz, Press & Schechter pro- 
posed the following analytic multiplicity function: 



(2) 



PS: uf(v) =W-i/exp(- 1 //2) 



This PS multiplicity function has successfully predicted the 
num erical m ultiplicity functi on in a qualitative manner. How- 
ever. lSheth & Tormenl JT9 99) proposed an alternative analytic 
multiplicity function that could better reproduce the numeri- 
cal multiplicity function: 

ST: vf(v)=A(\ + v'- 2p )^v'^(-v' 2 /2), (3) 



where v' 2 = av 2 , and A is a normalization factor defined so 
that L f(v)dv = 1 . Because of this unity constraint, A is not 
an independent parameter, but is expressed in forms of p: 



l+2- p n- 1/2 T(l/2-p) 



-1 



(4) 



ISheth & Tormlnl (l 999 ) gave the best-fit parameter values of 
a = 0.707,p = 0.3,andA = 0.322. The PS multiplicity function 
is included in this ST multiplicity function with a = l,p = 0, 
andA = l/2. 

The ST multiplicity function has a maximum at v = ^max 
that satisfies the following equation: 



2p 



max 



+ \)(v' r 



max 



D + 2/?=0, 



(5) 



where ^ ma x 2 = flt/ max 2 - It is trivial to see z^max = 1 for the 
c ase o f the PS multiplicity function. 

J01 also proposed an analytic multiplicity function which 
gives a fit to their numerical multiplicity function: 

J01: ^/(>) = 0.315exp(-|0.61+logzy-log6 e | 3 ' 8 ). (6) 

The valid range of this fitting formula is -1 .2 < log i^-log 5 C < 
1.05. Hereafter we assume that S c in the above equation is 
constant and takes the value in the Einstein-de Sitter universe, 
i.e. 5 C = (3/20)(12tt) 2 / 3 - 1.686 in order to compare the J01 
multiplicity function with other functions. 



3. SIMULATIONS AND RESULTS 

We used the Adaptive Mesh Refinement A^-body code de- 
veloped by Yahagi (2002), which is a ve ctorized and par- 
alleliz ed version of the code described in Yahagi & Yoshiil 
(2001). All five runs of simulations we performed adopt 
the ACDM cosmological parameters of il m = 0.3, fl\ = 0.7, 
h = 0.7, and us = 1.0, using 512 3 particles in common. The 
size of the finest mesh is 1/64 of the base mesh, and the force 
dynamic range is 2 15 = 32768. Other simulation parameters, 
such as the box size and the particle mass are given in Tabled 
Initial conditio ns were gener ated by the GRAFIC2 code pro- 
vided by Bertschinger ( 2001 ) using the power spectrum given 
by Bardeen et al. ( 1986). Five runs produced 529 files, and 
each of them was analyzed by the FoF algorithm with a con- 
stant linking length of 0.2 times the mean particle separation. 
Details of the simulations will be given in Yahagi et al., in 
preparation. 

The mass functions of dark halos are shown in Figure ^ 
The upper and lower panels show the mass functions at z = 3 
and z = 0, respectively. The PS mass function is represented 
by solid lines, and the ST mass function by dashed lines. At 
both redshifts, the numerical mass functions from our simu- 
lations agree with the ST mass function in a mass range of 
10 10 M Q <M< 10 I3 M Q . 

The numerical multiplicity functions are shown by crosses 
in five panels of Figure|2] All the data from the initial redshift 
to the present z = is compiled to draw the average curves 
(crosses) with error bars indicating the epoch to epoch varia- 
tion. In the panel (f), all the numerical multiplicity functions 
are shown by thin lines. Dark halos that consist of less than 
600 particles are not used in calculating the multiplicity func- 
tion, and 1/64 dex-sized bins containing less than 100 halos 
are excluded to avoid the contamination of the rare objects. 
Three analytic multiplicity functions described in the previ- 
ous section are also shown in this figure, that is PS (solid 
lines), ST (dashed lines), and J01 (dotted lines). The best- 
fit functions based on the ST functional form are also shown 
by dot-dashed lines. The best-fit parameters are given in Ta- 
ble 13 and are very close to those of W02> Since the data are 
available only in the region at v < 3, these functions could be 
erroneous at v > 3. 

In most cases, the numerical multiplicity functions and the 
best-fit functions to them are consistent with the ST and J01 
multiplicity functions at v > 3. However, each of the numeri- 
cal multiplicity functions reside between the ST and J01 func- 
tions at 1.5 < v < 3, and is below the ST function at v < 1 
except for the 35b run. The numerical multiplicity functions 
have an apparent peak at v ~ 1, instead of a plateau as seen in 
the J01 function. We here proposed the following function to 
fit to the numerical multiplicity function: 

vf(v) =A[l+(Bv/V2) c ]is D exp[-(Bv/V2) 2 ], (7) 

where, A is a normalization factor to satisfy the unity con- 
straint, J. f(v)dv = 1, therefore 



A = 2(B/V2) D {T[D/2] + r[(C+D)/2]Y 



(8) 



The best-fit parameters are given as B=0.893, C=1.39, and 
D=0.408, and from these parameters, A is constrained so that 
A = 0.298. This best-fit function from equationQis shown in 
Figure[3]and is only valid at 0.3 < v < 3. 

Since the unity constraint is only an assumption, we can 
relax this constraint and treat A as a free parameter. Actu- 
ally, almost a half of the particles are not bound to any ha- 
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los as shown in Table [3] On the other hand, the best-fit pa- 
rameters without the unity constraint using equation [3] are 
A = 0.320, B = 0.664, C = 1.99, and D = 0.36, for which 
!™ f(v)dv =1.2713. Since the fraction of bound particle 
should not exceed unity, we assign the unity constraint to all 
the best-fit functions below. 

We also investigate the peak position of the multiplicity 
functions. Since deriving fmax directly from the numer- 
ical multiplicity function is difficult due to the numerical 
fluctuation around the maximum value, we derived ^max 
from the best fit ST function by solving equation [5] These 
t'max are given in Table 12 and are very close to unity. For 
reference, fmax = 0.916 for the J01 multiplicity function, 
and t^nax = 0.881 for th e multiplicity function proposed by 
lLee & Shandarhl (fT998Tl . 

We also checked the time dependence of the multiplicity 
function. Figure |4] shows the multiplicity function from the 
35a run, for four redshift ranges of < z < 1 (circles), 1 < 
z < 3 (squares), 3 < z < 6 (triangles), and z > 6, (crosses). 
At high redshifts, high-^ halos in the exponential part of the 
best-fit ST function and equation are probed. As redshift 
decreases, the probe window moves to the lower-f region. 

Since there are no modes of density fluctuations whose 
wave length is larger than the box size on which the periodic 
boundary is placed. The numerical multiplicity function is 
conditional such that f(v\5box = 0), where 8b x is the density 
contrast smoothed over a mass scale comparable to the box 
size. Since our box size is smaller than that of other groups, 
we have estimated this box-size effect comparing the uncon- 
ditional PS multiplicity functi on using the sharp fc-space fil- 
ter with the con ditional one (Bond et ali i 1 99 It iBowerlll 99 It 
lLacev & ColeTl 993). 



take values between the ST and J01 functions. Although these 
differences are within 1 a error bars, they are pos sibly due to 

i-rr . , . . T77, r; "IjinnM 



Vfc 



2 v-a>bo. 



7T(T 



.£2)3/2 



exp 



(v-€V box ) 2 

' 2(1 -e 2 ) 



(9) 



where e = (Jboxl '■ Setting Vb ox = 0, this effect is found to be not 
so large at v < 3 for the box size of 35/i _1 Mpc, and it makes 
the universality of the multiplicity function even worse. 

4. DISCUSSION 

We have performed five runs of A^-body simulations with 
high mass resolution in order to study the behavior of the nu- 
merical multiplicity function in the low-mass range, or the 
low-i/ region. Throughout the peak range of, 0.3 < v < 3, the 
ST functional form provides a good fit to them with parameter 
values of a = 0.664, p = 0. 321, a nd A = 0.301. These values 
are very close to those of W02. Our numerical multiplicity 
functions have a peak at v ~ 1 as in the ST function, instead 
of a plateau in the J01 function. 

However, some detailed discrepancies are seen between our 
numerical multiplicity functions and the ST and J01 analytic 
functions. First, in the low-^ region of v < 1, our numeri- 
cal multiplicity functions systematically fall below the ST and 
the J01 functions, while they are consistent with that of lW02l 
On the other hand, in the high-^ region, where v is signifi- 
cantly larger than unity, our numerical multiplicity functions 



the different box sizes adopted. Sheth & Tormen ( 1999) used 
the data from the GIF simulations (Kauffman et al. 1999) with 
the box size of 144/r'Mpc or less, and W02 mainly used the 
data from the box size of 200/i _1 Mpc. On the other hand, 
iJOll used the 3000 /i _I Mpc simulation at maximum. We have 
taken into account the box-size effect using the conditional PS 
multiplicity function (equation |9j but this effect is found to 
be too small to resolve this problem. Introducing the estima- 
tion using the conditional multiplicity function based on the 
unconditional multiplicity function which fits the numerical 
multiplicity function well, such as the ST multiplicity func- 
tion, would resolve this problem. However, such an improved 
estimation of the box-size effect might be weaker than the es- 
timation based on the PS function, because the unconditional 
ST function at v ~ 1 has a broad peak that is below that of the 
PS function. The fact that our numerical multiplicity func- 
tions keep the universality supports this line of argument. 

Thus, there are two discrepancies remained. One is the 
discrepancy between the numerical and analytical multiplic- 
ity functions. Although our newly proposed functional form 
(equation[7} provides a better fit when compared with the ST 
functional form (equation [3}, we need an analytic function 
based on a theoretical background which fits the numerical 
multiplicity function even better. The other one is the dis- 
crepancy in the numerical multiplicity functions from various 
simulation runs. There are three strategies to resolve this dis- 
crepancy. The first is to run simulations having still higher 
mass dynamic range free from the box size effect . The second 
is to increase the number of realizations as W02 did, because 
there is a scatter from the runs using the same box size. The 
third is to run simulations whose box size is smaller than that 
of the present work, although it might sound contradictorily. 
From simulations with smaller box size, we will obtain the in- 
formation on the conditional multiplicity function which co- 
incides with the unconditional multiplicity function at v <C 1 . 
Comparing the unconditional multiplicity function from sim- 
ulations with a large box size and the conditional multiplicity 
function from those of a small box size will offer not only 
the clues to resolve the above mentioned discrepancies, but 
also insights into the mechanism how the PS ansatz works to 
reproduce the numerical multiplicity function. 
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FIG. 1 . — The mass functions of dark halos at (a) z = 3 and (b) z = 0. The ST mass function (dashed line) agrees with the numerical mass functions (symbols) 
fairly well at both redshifts, while the PS mass function (solid line) does not agree with them. 
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FIG. 2. — The numerical multiplicity functions from five runs of our simulations (crosses with error bars) are shown in the panels (a-e), except for the panel 
(/) which shows the results from all the runs (thin lines). The numerical multiplicity functions are derived by compiling all the data from the initial redshift to 
the present z = 0. Crosses and error bars represent their average and rms, respectively. Also shown in each panel are the PS multiplicity function (solid line), 
the ST multiplicity function (dashed line), the J01 multiplicity function (dotted line), and the best-fit function using the ST functional form (dot-dashed line) for 
which adopted parameters are given in Tablel2l The J01 multiplicity function, originally given as a function of a i Jenkins et al. 2001 ), is expressed here in terms 
of v = &c/cr, assuming that 8 C is constant although it varies slightly in the ACDM universe. In the high-mass range (v >1), the numerical multiplicity functions 
reside between the ST and J01 functions, and its maximum value at v ~ 1 is below those of the ST and J01 functions. 



Mass function of low mass dark halos 



7 



0.35 




V 

FIG. 3. — The best-fit multiplicity function using equation Q(so/;'d line). This function well represents the numerical multiplicity functions indicated by 
symbols with error bars. Also shown for comparison are the ST (dashed line), J01 (dotted line), and best-fit ST (dot-dashed line) functions. For 1.5 < v < 3, all 
the numerical functions reside between the ST and the J01 functions. Even the best-fit ST function cannot represent the multiplicity functions well in this region. 
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FIG. 4. — Time dependence of the multiplicity function from the 35a run; < z < 1 (circles), 1 < z < 3 (squares), 3 < z < 6 (triangles), and z > 6 (crosses). 
Also shown are the best fit function for all the runs using the ST functional form (dot-dashed line) and equationQ(so/rcf line). At high redshifts, high-i/ halos in 
the exponential part of the ST and equationQare probed. As redshift decreases, the probe window moves to the lower-!/ region. 
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Table 1. Simulation Parameters 



Model 


L [/r'Mpc] 


m ptc i [h 'M ] 


Zstart 


35a.. 


35 


2.67 X 10 7 


50 


35b.. 


35 


2.67 X 10 7 


50 


70a . . 


70 


2.13 X 10 s 


41 


70b . . 


70 


2.13 X 10 s 


41 


140. . 


140 


1.70 X 10 9 


33 



NOTE. — All the simulations adopt the cosmological parameters of fl„ = 0.3, O a = 0.7, h = 0.7, and <x 8 = 1.0. The number of particles is 512 3 and the force dynamic range of 
simulations is 2 15 = 32768 in common. L is the box size, m ptd is the particle mass, and z sta rt i s me initial redshift. 
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Table 2. Best-Fit Parameters 



Model 


CI 


p 


A 


fmax 


35a . . 


0.665 


0.317 


0.305 


0.998 


35b.. 


0.715 


0.303 


0.320 


0.975 


70a.. 


0.666 


0.327 


0.293 


0.988 


70b.. 


0.614 


0.325 


0.296 


1.031 


140.. 


0.658 


0.321 


0.300 


0.999 


all . . . 


0.664 


0.321 


0.301 


0.996 


PS... 


1 





0.5 


1 


ST... 


0.707 


0.3 


0.322 


0.983 


W02. 


0.64 


0.34 




0.995 



NOTE. — The ST functional form in equation's used to fit to the numerical multiplicity function. 
a The value of v at which the best-fit ST function attains a maximum. 
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Table 3. Fraction of Unbound Particles at 



Model 


*ub" 


35a.. 


0.472 


35b.. 


0.472 


70a.. 


0.530 


70b . . 


0.528 


140.. 


0.601 



"Fraction of unbound particles. 



